Flowchart

library(GENIE3)
library(doParallel)
library(igraph)
library(tidyverse)
library(DT)
library(reticulate)
library(learn2count)
library(rbenchmark)
library(reshape2)
library(gridExtra)
library(DiagrammeR)
library(pROC)
library(JRF)
library(DiagrammeRsvg)
library(rsvg)
library(RColorBrewer)
library(rbenchmark)
library(ZILGM)
library(patchwork)
library(scales)
library(INetTool)

use_python("/usr/bin/python3", required = TRUE)
arboreto <- import("arboreto.algo")
pandas <- import("pandas")
numpy <- import("numpy")

execution_times <- list()

source("generate_adjacency.R")
source("symmetrize.R")
source("pscores.R")
source("plotg.R")
source("compare_consensus.R")
source("create_consensus.R")
source("earlyj.R")
source("plotROC.R")
source("cutoff_adjacency.R")
source("infer_networks.R")

grViz_output <- DiagrammeR::grViz("
digraph biological_workflow {
  # Set up the graph attributes
  graph [layout = dot, rankdir = TB]

  # Define consistent node styles
  node [shape = rectangle, style = filled, color = lightblue, fontsize = 12]

  # Define nodes for each step
  StartNode [label = 'Ground Thruth - String Regulatory Network', shape = oval, color = seagreen, fontcolor = black]
  AdjacencyMatrix [label = 'Thruth Adjacency Matrix', shape = rectangle, color = seagreen]
  SimulateData [label = 'Simulate Single-Cell Data', shape = rectangle, color = goldenrod]

  # Reconstruction using Three Packages
  LateIntegration [label = 'Late\nIntegration', shape = oval, color = khaki]
  EarlyIntegration [label = 'Early\nIntegration', shape = oval, color = khaki]
  Jointanalysis [label = 'Joint\nanalysis', shape = oval, color = khaki]
  

  # Process 
  earlyj [label = 'earlyj.R', shape=diamond, color=lightblue, fontcolor=black]
  networkinference [label = 'infer_networks.R\nGENIE3\nGRNBoost2\nZILGM\nJRF', shape = rectangle, color = goldenrod, fontcolor=black]
  symmetrize [label = 'symmetrize.R', shape = rectangle, color = goldenrod, fontcolor=black]
  plotROC [label = 'plotROC.R', shape=diamond, color=lightblue, fontcolor=black]
  generateadjacency [label='generate_adjacency.R\nWeighted Adjacency', shape=rectangle, color=goldenrod, fontcolor=black]
  cutoffadjacency [label='cutoff_adjacency.R\nBinary Adjacency', shape=rectangle, color=goldenrod, fontcolor=black]
  pscores [label='pscores.R\nTPR\nFPR\nF1\nAccuracy\nPrecision', shape=diamond, color=lightblue, fontcolor=black]
  voting [label='Voting\nUnion\nIntersection', shape=diamond, color=lightblue, fontcolor=black]
  plotgcompare  [label='plotg.R\ncompare_consesus.R\nPlot Graphs', shape=rectangle, color=goldenrod, fontcolor=black]

  # Define the workflow structure
  StartNode -> AdjacencyMatrix
  AdjacencyMatrix -> SimulateData
  SimulateData -> LateIntegration
  SimulateData -> EarlyIntegration
  SimulateData -> Jointanalysis
  EarlyIntegration -> earlyj
  earlyj -> networkinference
  LateIntegration -> networkinference
  Jointanalysis -> networkinference
  networkinference -> symmetrize
  symmetrize -> plotROC
  symmetrize -> generateadjacency
  generateadjacency -> cutoffadjacency
  cutoffadjacency -> pscores
  cutoffadjacency -> voting
  voting -> plotgcompare
  voting -> pscores
}
")

svg_code <- export_svg(grViz_output)
rsvg::rsvg_png(charToRaw(svg_code), "./../analysis/flowchart.png")

grViz_output

Tcell Ground Truth

adjm <- read.table("./../data/adjacency_matrix.csv", header = T, row.names = 1, sep = ",") %>% as.matrix()
diag(adjm) <- 0

adjm %>%
    datatable(extensions = 'Buttons',
            options = list(
              dom = 'Bfrtip',
              buttons = c('csv', 'excel'),
              scrollX = TRUE,
              pageLength = 10), 
            caption = "Ground Truth")

gtruth <- igraph::graph_from_adjacency_matrix(adjm, mode = "undirected", diag = F)

num_nodes <- vcount(gtruth)
num_edges <- ecount(gtruth)

set.seed(1234)
plot(gtruth, 
     main = paste("Ground Truth\nNodes:", num_nodes, "Edges:", num_edges),
     vertex.label.color = "black",
     vertex.size = 6, 
     edge.width = 2, 
     vertex.label = NA,
     vertex.color = "steelblue",
     layout = igraph::layout_with_fr)

Simulate Data

ncell <- 500
nodes <- nrow(adjm)

set.seed(1130)
mu_values <- c(3, 6, 9)
theta_values <- c(1, 0.7, 0.5)

count_matrices <- lapply(1:3, function(i) {
  set.seed(1130 + i)
  mu_i <- mu_values[i]
  theta_i <- theta_values[i]
  
  count_matrix_i <- simdata(n = ncell, p = nodes, B = adjm, family = "ZINB", 
                            mu = mu_i, mu_noise = 1, theta = theta_i, pi = 0.2)
  
  count_matrix_df <- as.data.frame(count_matrix_i)
  colnames(count_matrix_df) <- colnames(adjm)
  rownames(count_matrix_df) <- paste("cell", 1:nrow(count_matrix_df), sep = "")
  
  return(count_matrix_df)
})

count_matrices[[1]] %>%
    datatable(extensions = 'Buttons',
            options = list(
              dom = 'Bfrtip',
              buttons = c('csv', 'excel'),
              scrollX = TRUE,
              pageLength = 10), 
            caption = "Simulated count matrix")

saveRDS(count_matrices, "./../analysis/count_matrices.RDS")

Matrices Integration

Late Integration

GENIE3

set.seed(1234)
tictoc::tic("GENIE3 late 15Cores")
genie3_late <- infer_networks(count_matrices, method="GENIE3")
saveRDS(genie3_late, "./../analysis/genie3_late.RDS")
execution_times[['GENIE3 late 15Cores']] <- tictoc::toc(log = TRUE)$toc[[1]]

## GENIE3 late 15Cores: 122.844 sec elapsed

genie3_late[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 output")

symmetrize Output and ROC

genie3_late_wadj <- generate_adjacency(genie3_late, ground.truth = adjm)
sgenie3_late_wadj <- symmetrize(genie3_late_wadj, weight_function = "mean")
genie3_late_auc <- plotROC(sgenie3_late_wadj, adjm, plot_title = "ROC curve - GENIE3 Late Integration")

sgenie3_late_wadj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 symmetrize output")

Generate Adjacency and Apply Cutoff

sgenie3_late_adj <- cutoff_adjacency(count_matrices = count_matrices,
                 weighted_adjm_list = sgenie3_late_wadj, 
                 ground.truth = adjm,
                 n = 3,
                 method = "GENIE3")

## Matrix 1 Mean 95th Percentile Cutoff: 0.009591004 
## Matrix 2 Mean 95th Percentile Cutoff: 0.009617292 
## Matrix 3 Mean 95th Percentile Cutoff: 0.009778179

sgenie3_late_adj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 adjacency")

Comparison with the Ground Truth

scores.genie3.late.all <- pscores(adjm, sgenie3_late_adj)

scores.genie3.late.all$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

plots <- plotg(sgenie3_late_adj)

consesusm <- create_consensus(sgenie3_late_adj, method="vote")
consesusu <- create_consensus(sgenie3_late_adj, method="union")
consesunet <- create_consensus(adj_matrix_list = sgenie3_late_adj, weighted_list = sgenie3_late_wadj, method = "INet", threshold = 0.05, ncores = 15)

## [1] 0.2726247
## [1] 0.06366566

scores.genie3.late <- pscores(adjm, list(consesusm))

scoresu.genie3.late <- pscores(adjm, list(consesusu))

scoresnet.genie3.late <- pscores(adjm, list(consesunet))

scores.genie3.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores vote")

scores.genie3.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores union")

scoresnet.genie3.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores INet")

ajm_compared <- compare_consensus(consesusm, adjm)

ajm_compared <- compare_consensus(consesusu, adjm)

ajm_compared <- compare_consensus(consesunet, adjm)

GRNBoost2

set.seed(1234)
tictoc::tic("GRNBoost2 late")
grnb_late <- infer_networks(count_matrices, method="GRNBoost2")
saveRDS(grnb_late, "./../analysis/grnb_late.RDS")
execution_times[['GRNBoost2 late']] <- tictoc::toc(log = TRUE)$toc[[1]]

## GRNBoost2 late: 8.639 sec elapsed

grnb_late[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 output")

symmetrize Output and ROC

grnb_late_wadj <- generate_adjacency(grnb_late, ground.truth = adjm)
sgrnb_late_wadj <- symmetrize(grnb_late_wadj, weight_function = "mean")
grnb_late_auc <- plotROC(sgrnb_late_wadj, adjm, plot_title = "ROC curve - GRNBoost2 Late Integration")

sgrnb_late_wadj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 symmetrize output")

Generate Adjacency and Apply Cutoff

sgrnb_late_adj <- cutoff_adjacency(count_matrices = count_matrices,
                 weighted_adjm_list = sgrnb_late_wadj, 
                 ground.truth = adjm,
                 n = 3,
                 method = "GRNBoost2")

## Matrix 1 Mean 95th Percentile Cutoff: 0.8642676 
## Matrix 2 Mean 95th Percentile Cutoff: 0.845515 
## Matrix 3 Mean 95th Percentile Cutoff: 0.8459939

sgrnb_late_adj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 adjacency")

Comparison with the Ground Truth

scores.grn.late.all <- pscores(adjm, sgrnb_late_adj)

scores.grn.late.all$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

plots <- plotg(sgrnb_late_adj)

consesusm <- create_consensus(sgrnb_late_adj, method="vote")
consesusu <- create_consensus(sgrnb_late_adj, method="union")
consesunet <- create_consensus(adj_matrix_list = sgrnb_late_adj, weighted_list = sgrnb_late_wadj, method = "INet", threshold = 0.05, ncores = 15)

## [1] 0.2747324
## [1] 0.03408968

scores.grn.late <- pscores(adjm, list(consesusm))

scoresu.grn.late <- pscores(adjm, list(consesusu))

scoresnet.grn.late <- pscores(adjm, list(consesunet))

scores.grn.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores vote")

scoresu.grn.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores union")

scoresnet.grn.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores INet")

ajm_compared <- compare_consensus(consesusm, adjm)

ajm_compared <- compare_consensus(consesusu, adjm)

ajm_compared <- compare_consensus(consesunet, adjm)

ZILGM Park

set.seed(1234)
tictoc::tic("ZILGM late 15Cores")
zilgm_late <- infer_networks(count_matrices_list = count_matrices, method = "ZILGM", adjm = adjm)

## learning for NBII graphical model 
## nlambda : 50
## penalty function : LASSO
## update type : IRLS
## Conducting sampling in progress :  10 % 
Conducting sampling in progress :  20 % 
Conducting sampling in progress :  30 % 
Conducting sampling in progress :  40 % 
Conducting sampling in progress :  50 % 
Conducting sampling in progress :  60 % 
Conducting sampling in progress :  70 % 
Conducting sampling in progress :  80 % 
Conducting sampling in progress :  90 % 
Conducting sampling in progress :  100 % 
learning for NBII graphical model 
## nlambda : 50
## penalty function : LASSO
## update type : IRLS
## Conducting sampling in progress :  10 % 
Conducting sampling in progress :  20 % 
Conducting sampling in progress :  30 % 
Conducting sampling in progress :  40 % 
Conducting sampling in progress :  50 % 
Conducting sampling in progress :  60 % 
Conducting sampling in progress :  70 % 
Conducting sampling in progress :  80 % 
Conducting sampling in progress :  90 % 
Conducting sampling in progress :  100 % 
learning for NBII graphical model 
## nlambda : 50
## penalty function : LASSO
## update type : IRLS
## Conducting sampling in progress :  10 % 
Conducting sampling in progress :  20 % 
Conducting sampling in progress :  30 % 
Conducting sampling in progress :  40 % 
Conducting sampling in progress :  50 % 
Conducting sampling in progress :  60 % 
Conducting sampling in progress :  70 % 
Conducting sampling in progress :  80 % 
Conducting sampling in progress :  90 % 
Conducting sampling in progress :  100 % 
Assigned lamb to lambda_results[[ 1 ]] with 50 matrices.
## Assigned lamb to lambda_results[[ 2 ]] with 50 matrices.
## Assigned lamb to lambda_results[[ 3 ]] with 50 matrices.

saveRDS(zilgm_late, "./../analysis/zilgm_late.RDS")
execution_times[['ZILGM late 15Cores']] <- tictoc::toc(log = TRUE)$toc[[1]]

## ZILGM late 15Cores: 3349.755 sec elapsed

est_graphs <- zilgm_late$network_results
lambdas <- zilgm_late$lambda_results

plotROC(lambdas[[1]], adjm, plot_title = "Lambda ROC Matrix1", is_binary = T)

plotROC(lambdas[[2]], adjm, plot_title = "Lambda ROC Matrix2", is_binary = T)

plotROC(lambdas[[3]], adjm, plot_title = "Lambda ROC Matrix3", is_binary = T)

consensus_matrices <- vector("list", 50)

for (i in 1:50) {
  ranklambda <- list(lambdas[[1]][[i]], lambdas[[2]][[i]], lambdas[[3]][[i]])
  consensus_matrices[[i]] <- create_consensus(ranklambda, method="vote")
}

zilgm_late_auc <- plotROC(consensus_matrices, adjm, plot_title = "ROC voting lambda", is_binary = T)

Comparison with the Ground Truth

scores.zilgm.late.all <- pscores(adjm, est_graphs)

scores.zilgm.late.all$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

plots <- plotg(est_graphs)

consesusm <- create_consensus(est_graphs, method="vote")
consesusu <- create_consensus(est_graphs, method="union")

scores.zilgm.late <- pscores(adjm, list(consesusm))

scoresu.zilgm.late <- pscores(adjm, list(consesusu))

scores.zilgm.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

scoresu.zilgm.late$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

ajm_compared <- compare_consensus(consesusm, adjm)

ajm_compared <- compare_consensus(consesusu, adjm)

Early Integration

early_matrix <- list(earlyj(count_matrices))

GENIE3

set.seed(1234)
tictoc::tic("GENIE3 early 15Cores")
genie3_early <- infer_networks(early_matrix, method="GENIE3")
execution_times[['GENIE3 early 15Cores']] <- tictoc::toc(log = TRUE)$toc[[1]]

## GENIE3 early 15Cores: 162.813 sec elapsed

saveRDS(genie3_early, "./../analysis/genie3_early.RDS")

genie3_early[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 output")

symmetrize Output and ROC

genie3_early_wadj <- generate_adjacency(genie3_early, ground.truth = adjm)
sgenie3_early_wadj <- symmetrize(genie3_early_wadj, weight_function = "mean")
genie3_early_auc <- plotROC(sgenie3_early_wadj, adjm, plot_title = "ROC curve - GENIE3 Early Integration")

sgenie3_early_wadj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 symmetrize output")

Generate Adjacency and Apply Cutoff

sgenie3_early_adj <- cutoff_adjacency(count_matrices = early_matrix,
                 weighted_adjm_list = sgenie3_early_wadj, 
                 ground.truth = adjm,
                 n = 2,
                 method = "GENIE3")

## Matrix 1 Mean 95th Percentile Cutoff: 0.01016292

sgenie3_early_adj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GENIE3 adjacency")

Comparison with the Ground Truth

scores.genie3.early <- pscores(adjm, sgenie3_early_adj)

scores.genie3.early$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

plots <- plotg(sgenie3_early_adj)

ajm_compared <- compare_consensus(sgenie3_early_adj[[1]], adjm)

GRNBoost2

set.seed(1234)
tictoc::tic("GRNBoost2 early")
grnb_early <- infer_networks(early_matrix, method="GRNBoost2")
execution_times[['GRNBoost2 early']] <- tictoc::toc(log = TRUE)$toc[[1]]

## GRNBoost2 early: 12.767 sec elapsed

saveRDS(grnb_early, "./../analysis/grnb_early.RDS")

grnb_early[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 output")

symmetrize Output and ROC

grnb_early_wadj <- generate_adjacency(grnb_early, ground.truth = adjm)
sgrnb_early_wadj <- symmetrize(grnb_early_wadj, weight_function = "mean")
grnb_early_auc <- plotROC(sgrnb_early_wadj, adjm, plot_title = "ROC curve - GRNBoost2 Early Integration")

grnb_early_wadj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 symmetrize output")

Generate Adjacency and Apply Cutoff

sgrnb_early_adj <- cutoff_adjacency(count_matrices = early_matrix,
                 weighted_adjm_list = sgrnb_early_wadj, 
                 ground.truth = adjm,
                 n = 2,
                 method = "GRNBoost2")

## Matrix 1 Mean 95th Percentile Cutoff: 4.228947

sgrnb_early_adj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "GRNBoost2 adjacency")

Comparison with the Ground Truth

scores.grn.early <- pscores(adjm, sgrnb_early_adj)

scores.grn.early$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

plots <- plotg(sgrnb_early_adj)

ajm_compared <- compare_consensus(sgrnb_early_adj[[1]], adjm)

ZILGM Park

set.seed(1234)
tictoc::tic("ZILGM early 15Cores")
zilgm_late <- infer_networks(count_matrices_list = early_matrix, method = "ZILGM", adjm = adjm)

## learning for NBII graphical model 
## nlambda : 50
## penalty function : LASSO
## update type : IRLS
## Conducting sampling in progress :  10 % 
Conducting sampling in progress :  20 % 
Conducting sampling in progress :  30 % 
Conducting sampling in progress :  40 % 
Conducting sampling in progress :  50 % 
Conducting sampling in progress :  60 % 
Conducting sampling in progress :  70 % 
Conducting sampling in progress :  80 % 
Conducting sampling in progress :  90 % 
Conducting sampling in progress :  100 % 
Assigned lamb to lambda_results[[ 1 ]] with 50 matrices.

saveRDS(zilgm_late, "./../analysis/zilgm_early.RDS")
execution_times[['ZILGM early 15Cores']] <- tictoc::toc(log = TRUE)$toc[[1]]

## ZILGM early 15Cores: 3122.533 sec elapsed

est_graphs <- zilgm_late$network_results
lambdas <- zilgm_late$lambda_results

zilgm_early_auc <- plotROC(lambdas[[1]], adjm, plot_title = "ROC ZILGM early", is_binary = T)

Comparison with the Ground Truth

scores.zilgm.early <- pscores(adjm, est_graphs)

scores.zilgm.early$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

ajm_compared <- compare_consensus(est_graphs[[1]], adjm)

Joint Integration

Joint Random Forest

#https://cran.r-project.org/src/contrib/Archive/JRF/
#install.packages("/home/francescoc/Downloads/JRF_0.1-4.tar.gz", repos = NULL, type = "source")
set.seed(1234)
tictoc::tic("JRF")
jrf_mat <- infer_networks(count_matrices, method="JRF")
execution_times[['JRF']] <- tictoc::toc(log = TRUE)$toc[[1]]

## JRF: 1171.366 sec elapsed

jrf_mat[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "JRF output")

Their Permutation

jrf_matrices <- lapply(count_matrices, t)
jrf_matrices_norm <- lapply(jrf_matrices,function(x) {
  (x - mean(x)) / sd(x)
  })

genes <- rownames(jrf_matrices_norm[[1]])

set.seed(1234)
out.perm <- Run_permutation(jrf_matrices_norm,mtry=round(sqrt(length(genes)-1)),ntree=500, genes,5)

myJRF_network <- function(out.jrf, out.perm, TH) {
  
  nclasses <- dim(out.perm)[3]
  M <- dim(out.perm)[2]
  out <- vector("list", nclasses)  

  for (net in 1:nclasses) { 
    j.np <- sort(out.jrf[, 2 + net], decreasing = TRUE)
    FDR <- matrix(0, dim(out.perm)[1], 1)
    
    th <- NULL  
    for (s in 1:length(j.np)) { 
      FP <- sum(sum(out.perm[, , net] >= j.np[s])) / M
      FDR[s] <- FP / s
      
      if (FDR[s] > TH) {
        th <- j.np[s]
        break
      }
    }
    
    out[[net]] <- out.jrf[out.jrf[, 2 + net] > th, seq(1, 2)]
  }
  
  return(out)
}

mynet <- myJRF_network(jrf_mat[[1]],out.perm,0.05)
mynet

## [[1]]
##      gene1 gene2
## 5444 PTPRC CXCR4
## 
## [[2]]
##         gene1  gene2
## 2321    CXCR4  CCR10
## 2733     FCMR   CCR7
## 3308     PRF1    CD2
## 3842    TPCN1   CD38
## 3856      CD6   CD3D
## 3911  LDLRAP1   CD3D
## 4151    KLRB1   CD3G
## 4167     NKG7   CD3G
## 4921   ZNF683   CD8A
## 5444    PTPRC  CXCR4
## 5461     STAM  CXCR4
## 7062   PIK3R5    FGR
## 7224    ITGB2   FLNA
## 7291    FOXP3  FOXP1
## 7731    USP15   GLUL
## 7748    IL2RB   GNLY
## 8008    KLRB1   GZMB
## 8082    IL2RB   GZMH
## 8118    PTPRC   GZMH
## 8438    RACK1  IGF1R
## 9251     TFRC  ITGB1
## 9333    KLRF1  KLRB1
## 9420     PRF1  KLRD1
## 9669    PTPRC LGALS1
## 10883 RASGRP2  RAP1B
## 11178   SNX20 SELPLG
## 
## [[3]]
## [1] gene1 gene2
## <0 rows> (or 0-length row.names)

jrf_adj <- function(df_list, adjm) {
  # Ensure the genes (row and column names) are in the adjacency matrix
  genes <- rownames(adjm)
  
  # Initialize a list to store adjacency matrices
  adjacency_matrices <- list()
  
  # Loop through each data frame in the list
  for (i in seq_along(df_list)) {
    df <- df_list[[i]]
    
    # Initialize a new adjacency matrix with zeros
    adj_matrix <- adjm * 0  # Reset to zero for each df
    
    # Update the adjacency matrix for the gene pairs
    for (j in seq_len(nrow(df))) {
      gene1 <- df$gene1[j]
      gene2 <- df$gene2[j]
      
      # Check if both genes are present in the adjacency matrix
      if (gene1 %in% genes && gene2 %in% genes) {
        idx1 <- which(genes == gene1)  # Find the row index for gene1
        idx2 <- which(genes == gene2)  # Find the column index for gene2
        
        # Update the adjacency matrix to set an edge
        adj_matrix[idx1, idx2] <- 1
        adj_matrix[idx2, idx1] <- 1  # Ensure symmetry
      } else {
        warning(paste("Gene pair not found in adjm:", gene1, "-", gene2))
      }
    }
    
    # Store the resulting adjacency matrix in the list
    adjacency_matrices[[i]] <- adj_matrix
  }
  
  # Return the list of adjacency matrices
  return(adjacency_matrices)
}


jrf_adj_mynet <- jrf_adj(mynet, adjm)

jrf_auc_perm <- plotROC(jrf_adj_mynet, adjm, plot_title = "ROC curve - JRF perm out", is_binary = T)

scores.jrf.perm <- pscores(adjm, jrf_adj_mynet)

scores.jrf.perm$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

plots <- plotg(jrf_adj_mynet)

consesusu <- create_consensus(jrf_adj_mynet, method="union")

scoresu.jrf.perm <- pscores(adjm, list(consesusu))

scoresu.jrf.perm$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

ajm_compared <- compare_consensus(consesusu, adjm)

Prepare the output

jrf_list <- list()

importance_columns <- grep("importance", names(jrf_mat[[1]]), value = TRUE)

for (i in seq_along(importance_columns)) {
  # Select the 'gene1', 'gene2', and the current 'importance' column
  df <- jrf_mat[[1]][, c("gene1", "gene2", importance_columns[i])]
  
  # Rename the importance column to its original name (e.g., importance1, importance2, etc.)
  names(df)[3] <- importance_columns[i]
  
  # Add the data frame to the output list
  jrf_list[[i]] <- df
}

saveRDS(jrf_list, "./../analysis/jrf.RDS")

jrf_list[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "JRF output")

symmetrize Output and ROC

jrf_wadj <- generate_adjacency(jrf_list, ground.truth = adjm)
sjrf_wadj <- symmetrize(jrf_wadj, weight_function = "mean")
jrf_auc_mine <- plotROC(sjrf_wadj, adjm, plot_title = "ROC curve - JRF Late Integration")

sjrf_wadj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "JRF symmetrize output")

Generate Adjacency and Apply Cutoff

sjrf_adj <- cutoff_adjacency(count_matrices = count_matrices,
                 weighted_adjm_list = sjrf_wadj, 
                 ground.truth = adjm,
                 n = 3,
                 method = "JRF")

## Matrix 1 Mean 95th Percentile Cutoff: 4.699322 
## Matrix 2 Mean 95th Percentile Cutoff: 4.74195 
## Matrix 3 Mean 95th Percentile Cutoff: 4.790277

sjrf_adj[[1]] %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "JRF adjacency")

Comparison with the Ground Truth

scores.jrf.all <- pscores(adjm, sjrf_adj)

scores.jrf.all$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores")

plots <- plotg(sjrf_adj)

consesusm <- create_consensus(sjrf_adj, method="vote")
consesusu <- create_consensus(sjrf_adj, method="union")
consesunet <- create_consensus(adj_matrix_list = sjrf_adj, weighted_list = sjrf_wadj, method = "INet", threshold = 0.1, ncores = 15)

## [1] 0.3228302
## [1] 0.09162004

scores.jrf <- pscores(adjm, list(consesusm))

scoresu.jrf <- pscores(adjm, list(consesusu))

scoresnet.jrf <- pscores(adjm, list(consesunet))

scores.jrf$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores vote")

scoresu.jrf$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores union")

scoresnet.jrf$Statistics %>%
    datatable(extensions = 'Buttons',
              options = list(
                dom = 'Bfrtip',
                buttons = c('csv', 'excel'),
                scrollX = TRUE,
                pageLength = 10), 
              caption = "scores INet")

ajm_compared <- compare_consensus(consesusm, adjm)

ajm_compared <- compare_consensus(consesusu, adjm)

ajm_compared <- compare_consensus(consesunet, adjm)

#tictoc::tic()
#est <- PCzinb(as.matrix(count_matrices[[1]]), method="zinb1", maxcard=2)
#tictoc::toc()

#graph_est <- graph_from_adjacency_matrix(est, mode="undirected")
#plot(graph_est)

Method Comparison

time_data <- data.frame(
  Method = names(execution_times),
  Time_in_Hours = unlist(execution_times) / 3600
)
time_data$Time_in_Minutes <- time_data$Time_in_Hours * 60
time_data <- time_data[order(time_data$Time_in_Hours), ]
time_data$Method <- factor(time_data$Method, levels = time_data$Method)

time_data <- time_data %>%
  mutate(Method_Group = case_when(
    grepl("GENIE3", Method) ~ "GENIE3",
    grepl("GRNBoost2", Method) ~ "GRNBoost2",
    grepl("ZILGM", Method) ~ "ZILGM",
    grepl("JRF", Method) ~ "JRF"
  ))

method_colors <- c("GENIE3" = "darkblue", "GRNBoost2" = "darkgreen", "ZILGM" = "orange", "JRF" = "red")

time_plot <- ggplot(time_data, aes(x = Method, y = Time_in_Hours, fill = Method_Group)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = sprintf("%.1f min", Time_in_Minutes)), vjust = -0.5) +
  labs(title = "Execution Time for Each Method", y = "Time (Hours)", x = NULL) +
  scale_fill_manual(values = method_colors) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none"
  )

print(time_plot)

mplot <- list()

for (k in c("TPR", "FPR", "Precision", "F1", "Accuracy")) {
  mplot[[k]] <- data.frame(
    GENIE3_late = scores.genie3.late$Statistics[[k]],
    GENIE3_early = scores.genie3.early$Statistics[[k]],
    GRNBoost2_late = scores.grn.late$Statistics[[k]],
    GRNBoost2_early = scores.grn.early$Statistics[[k]],
    ZILGM_late = scores.zilgm.late$Statistics[[k]],
    ZILGM_early = scores.zilgm.early$Statistics[[k]],
    JRF = scores.jrf$Statistics[[k]]
  )
}

mplot[["AUC"]] <- data.frame(
    GENIE3_late = mean(genie3_late_auc$AUC),
    GENIE3_early = genie3_early_auc$AUC,
    GRNBoost2_late = mean(grnb_late_auc$AUC),
    GRNBoost2_early = grnb_early_auc$AUC,
    ZILGM_late = zilgm_late_auc$AUC,
    ZILGM_early = zilgm_early_auc$AUC,
    JRF = jrf_auc_mine$AUC
)

plot_data <- bind_rows(lapply(names(mplot), function(metric) {
  data.frame(
    Metric = metric,
    Method = names(mplot[[metric]]),
    Value = as.numeric(mplot[[metric]][1, ])
  )
}))

plot_data <- plot_data %>%
  mutate(Method_Group = case_when(
    grepl("GENIE3", Method) ~ "GENIE3",
    grepl("GRNBoost2", Method) ~ "GRNBoost2",
    grepl("ZILGM", Method) ~ "ZILGM",
    grepl("JRF", Method) ~ "JRF"
  )) %>%
  mutate(Method = factor(Method, levels = c(
    "GENIE3_early", "GENIE3_late", 
    "GRNBoost2_early", "GRNBoost2_late", 
    "ZILGM_early", "ZILGM_late", 
    "JRF"  # Ensure JRF is last
  )))

plots <- lapply(unique(plot_data$Metric), function(metric) {
  show_x_text <- metric %in% c("Accuracy", "AUC")
  
  ggplot(plot_data %>% filter(Metric == metric), aes(x = Method, y = Value, fill = Method_Group)) +
    geom_bar(stat = "identity", position = position_dodge(width = 0.8)) +
    labs(title = metric, y = "Value", x = NULL) +
    scale_y_continuous(limits = c(0, 1)) +  
    scale_fill_manual(values = method_colors) +
    theme_minimal() +
    theme(
      axis.text.x = if (show_x_text) element_text(angle = 45, hjust = 1) else element_blank(),
      axis.title.x = element_blank(),
      legend.position = "none"  
    )
})

final_plot <- (plots[[1]] | plots[[2]]) /
              (plots[[3]] | plots[[4]]) /
              (plots[[5]] | plots[[6]]) +
  plot_layout(guides = "collect") & 
  theme(legend.position = "bottom")

print(final_plot)

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=it_IT.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=it_IT.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=it_IT.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=it_IT.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] doRNG_1.8.2        rngtools_1.5.2     INetTool_0.1.0     scales_1.2.0      
##  [5] patchwork_1.1.1    ZILGM_0.1.1        RColorBrewer_1.1-3 rsvg_2.6.1        
##  [9] DiagrammeRsvg_0.1  JRF_0.1-4          pROC_1.18.0        DiagrammeR_1.0.11 
## [13] gridExtra_2.3      reshape2_1.4.4     rbenchmark_1.0.0   learn2count_0.3.2 
## [17] reticulate_1.34.0  DT_0.22            forcats_0.5.1      stringr_1.4.0     
## [21] dplyr_1.0.9        purrr_0.3.4        readr_2.1.2        tidyr_1.2.0       
## [25] tibble_3.1.7       ggplot2_3.3.6      tidyverse_1.3.1    igraph_2.0.3      
## [29] doParallel_1.0.17  iterators_1.0.14   foreach_1.5.2      GENIE3_1.16.0     
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.0                backports_1.4.1            
##   [3] plyr_1.8.7                  splines_4.1.0              
##   [5] crosstalk_1.2.0             GenomeInfoDb_1.30.1        
##   [7] digest_0.6.29               htmltools_0.5.2            
##   [9] fansi_1.0.3                 magrittr_2.0.3             
##  [11] tzdb_0.3.0                  modelr_0.1.8               
##  [13] matrixStats_0.62.0          colorspace_2.0-3           
##  [15] rvest_1.0.2                 haven_2.5.0                
##  [17] xfun_0.30                   crayon_1.5.1               
##  [19] RCurl_1.98-1.6              jsonlite_1.8.0             
##  [21] graph_1.72.0                survival_3.3-1             
##  [23] glue_1.6.2                  gtable_0.3.0               
##  [25] zlibbioc_1.40.0             distributions3_0.2.2       
##  [27] XVector_0.34.0              DelayedArray_0.20.0        
##  [29] V8_6.0.0                    car_3.0-13                 
##  [31] shape_1.4.6                 SingleCellExperiment_1.16.0
##  [33] BiocGenerics_0.40.0         abind_1.4-5                
##  [35] pscl_1.5.9                  DBI_1.1.2                  
##  [37] rstatix_0.7.0               Rcpp_1.0.8.3               
##  [39] datastructures_0.2.9        stats4_4.1.0               
##  [41] glmnet_4.1-8                htmlwidgets_1.5.4          
##  [43] httr_1.4.3                  WeightSVM_1.7-16           
##  [45] ellipsis_0.3.2              farver_2.1.0               
##  [47] pkgconfig_2.0.3             sass_0.4.1                 
##  [49] dbplyr_2.1.1                utf8_1.2.2                 
##  [51] labeling_0.4.2              tidyselect_1.1.2           
##  [53] rlang_1.1.4                 munsell_0.5.0              
##  [55] cellranger_1.1.0            tools_4.1.0                
##  [57] visNetwork_2.1.2            cli_3.3.0                  
##  [59] generics_0.1.2              statnet.common_4.10.0      
##  [61] broom_0.8.0                 evaluate_0.15              
##  [63] fastmap_1.1.0               yaml_2.3.5                 
##  [65] bst_0.3-24                  knitr_1.39                 
##  [67] fs_1.5.2                    caTools_1.18.2             
##  [69] tictoc_1.2.1                xml2_1.3.3                 
##  [71] mpath_0.4-2.26              multinet_4.2.1             
##  [73] compiler_4.1.0              rstudioapi_0.13            
##  [75] curl_4.3.2                  png_0.1-7                  
##  [77] ggsignif_0.6.3              reprex_2.0.1               
##  [79] bslib_0.3.1                 stringi_1.7.6              
##  [81] highr_0.9                   lattice_0.20-45            
##  [83] iZID_0.0.1                  Matrix_1.6-1.1             
##  [85] gbm_2.2.2                   vctrs_0.4.1                
##  [87] pillar_1.7.0                lifecycle_1.0.1            
##  [89] jquerylib_0.1.4             bitops_1.0-7               
##  [91] GenomicRanges_1.46.1        R6_2.5.1                   
##  [93] network_1.18.2              IRanges_2.28.0             
##  [95] codetools_0.2-18            MASS_7.3-57                
##  [97] assertthat_0.2.1            SummarizedExperiment_1.24.0
##  [99] withr_2.5.0                 S4Vectors_0.32.4           
## [101] GenomeInfoDbData_1.2.7      hms_1.1.1                  
## [103] grid_4.1.0                  rpart_4.1.16               
## [105] coda_0.19-4                 flux_0.3-0.1               
## [107] rmarkdown_2.14              MatrixGenerics_1.6.0       
## [109] carData_3.0-5               ggpubr_0.4.0               
## [111] numDeriv_2016.8-1.1         Biobase_2.54.0             
## [113] lubridate_1.8.0

scGRN testing

Francesco Cecere

2024-11-19

Flowchart

Tcell Ground Truth

Simulate Data

Matrices Integration

Late Integration

GENIE3

symmetrize Output and ROC

Generate Adjacency and Apply Cutoff

Comparison with the Ground Truth

GRNBoost2

symmetrize Output and ROC

Generate Adjacency and Apply Cutoff

Comparison with the Ground Truth

ZILGM Park

Comparison with the Ground Truth

Early Integration

GENIE3

symmetrize Output and ROC

Generate Adjacency and Apply Cutoff

Comparison with the Ground Truth

GRNBoost2

symmetrize Output and ROC

Generate Adjacency and Apply Cutoff

Comparison with the Ground Truth

ZILGM Park

Comparison with the Ground Truth

Joint Integration

Joint Random Forest

Their Permutation

Prepare the output

symmetrize Output and ROC

Generate Adjacency and Apply Cutoff

Comparison with the Ground Truth

Method Comparison